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PREFACE 


This  report  was  prepared  by  the  University  ol  Dayton  Research  Institute  (UDRI) 
for  the  US  Amy  Atmospheric  Sciences  Laboratory  (ASL)  during  the  period  1 
October  1980  to  31  March  1981.  The  report  describes  calculations  made  of  the 
aeroclynamlc  surface  roughness  for  a  region  of  Central  Europe.  The 
assumptions,  techniques,  and  mathematical  algorithms  used  for  this  task  are 
presented.  Two  sample  calculations  for  a  1  km  area  are  given. 

The  contract  technical  officer  for  this  work  was  John  T.  Marrs  of  ASL.  Work 
was  performed  at  UDRI  under  the  administrative  supervision  of  Nicholas  A. 
Engler.  The  technical  supervisor  was  James  K.  Luers.  Computer  programming 
and  data  processing  were  done  by  Nancy  J.  Fratlnl,  Jerry  6.  Jensen*  Steven  G. 
Vondrell,  and  Zalfa  A.  Abdelnour.  Pamela  S.  Ecker  edited  the  report  and 
Gretchen  Walther  prepared  the  typescript.  The  authors  Wish  to  acknowledge  the 
very  substantial  contributions  of  all  of  those  mentioned  above  to  the 
performance  of  this  project. 

This  particular  work  effort  was  suggested  by  Frank  V.  Hansen  who  also  supplied 
useful  Ideas  on  the  proper  approach  to  the  problem. 
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INTRODUCTION 


The  transport  and  diffusion  of  smoke  and  other  contaminants  In  the 
atmospheric  boundary  layer  defends  In  a  large  part  upon  the  aerodynamic  rough¬ 
ness  of  the  surface.  Surface  roughness  arises  from  the  small  scale  protru¬ 
sions  (roughness  elements)  which  penetrate  the  lowest  portions  of  the  boundary 
Tayer to-«xer£_a^trong  influence  on  the  mean  wind  profile  and  the  level  of 
turbulence.  Common  meteorological  practi ce  [ 1 ,2]  is  to  prescribe  the  roughness 
from  the  vegetation  (forests,  grasslands,  agriculture,  etc.),  the  man-made 
objects  (buildings,  cities),  or  the  other  surface  characteristics  (mud  flats, 
sea  surface,  etc.)  that  best  characterize  the  area  under  consideration.  The 
quantitative  measure  of  roughness  Is  the  roughness  length,  z  ,  which  Is  pro- 
portlonal  to.  If  not  a  direct  measure  of,  the  vertical  dimension  of  the 
Individual  roughness  elements.  Each  type  of  surface  may  be  assigned  a  value 
of  roughness  length— the  values  obtained  primarily  from  field  experiments. 

C 

Typical  values  range  from  zQ  *  10  m  for  mud  flats  or  Ice  to  zQ  “  10m  for  the 
tall  buildings  In  a  central  city. 

Field  studies  show  that  a  single  value  for  roughness  length  prescribed 
in  this  way  works  well  In  describing  the  mean  wind  and  turbulence  when  the 
terrain  Is  reasonably  flat  on  the  large  scale  and  of  a  homogeneous  character. 
However,  most  practical  situations  Involve  "complex  terrain"  where  the  eleva¬ 
tion  of  the  land  may  vary  considerably  over  relatively  short  distances  and  the 
surface  morphology  may  be  a  complicated  mixture  of  roughness  elements  that 
does  not  fall  Into  any  of  the  usual  tidy  classifications. 

Transport,  diffusion,  and  the  general  structure  of  the  boundary  layer  In 
regions  of  complex  terrain  Is  the  object  of  much  current  research.  Some 


[1]  Pasqulll,  F.f  Atmospheric  Diffusion,  (2nd  Edition),  Ellis 
Norwood  Ltd.,  i9/e. 

[2]  Hansen,  F.V.,  "Engineering  Estimates  for  the  Calculation  of 
Atmospheric  Dispersion  Coefficients,"  ASL  Internal  Report, 
U.S.  Army  Atmospheric  Sciences  Laboratory,  White  Sands  Army 
Missile  Range,  New  Mexico,  September  1979. 
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Investigators  [3,4]  have  suggested  that  the  general  elevation  changes  over  the 
scale  of  several  kilometers  or  more  my  contribute  to,  or  even  control,  the 
effective  surface  roughness  In  some  cases.  This  suggests  the  possibility  that 
roughness  lengths  can  be  computed  from  (sufficiently  detailed)  elevation  data 
alone  In  very  hilly  or  mountainous  terrain. 

This  report  describes  two  methods  of  computing  an  estimate  of  the  sur¬ 
face  roughness  from  detailed  terrain  descriptions  of  the  kind  that  may  be 
obtained  from  topographic  maps.  The  first  and  more  traditional  method,  tftlch 
we  call  the  Surface  Morphology  Roughness  (SMR)  estimation,  uses  descriptions 
of  vegetation  and  other  surface  features  to  obtain  roughness  lengths  without 
reference  to  the  elevation  field.  Since  seasonally  changing  plant  growth  Is 
often  the  major  surface  feature,  the  SMR  method  Includes  the  changes  In  zQ  by 
month  of  the  year.  The  second  method,  Terrain  Derived  Roughness  (TDR) 
estimation,  uses  just  the  elevation  data— assumed  to  be  available  on  a  square 
grid— to  compute  a  roughness  length  based  upon  changes  in  the  surface 
elevation.  Examples  of  each  computation  are  given  and  a  brief  comparison  of 
the  two  methods  Is  Included. 


[3]  Thompson,  R.  S.,  "Note  on  the  Aerodynamic  Roughness  Length 

for  Complex  Terrain,"  J.  Appl.  Meteor..  Vol  17.  1402-1403. 
(1978).  - 

[4]  Lettau,  H.,  "Note  on  Aerodynamic  Roughness-Parameter 
Estimation  on  the  Basis  of  Roughness -El ament  Description," 
J.  Appl.  Meteor..  Vol.  8,  828-832,  (1969). 
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STATEMENT  OF  THE  PROBLEM 

The  specific  problem  addressed  here  is  the  computation  of  roughness 
lengths  for  a  large  section  of  Central  Europe.  Two  sets  of  data  describing 
the  region  were  available.  The  first  and  smaller  set  consisted  of  elevation 
and  vegetation  data  supplied  on  a  100  meter  square  grid  covering  an  area  of  97 
by  125  kilometers.  For  brevity  we  call  this  area  region  A.  The  second  data 
set  consisted  of  only  elevation  data  but  covered  a  much  larger  area,  about 
400,000  km  .  We  call  this  region  B.  Data  specification  In  region  B  Is  also 
on  a  square  grid  but  with  a  spacing  of  63.5  meters  between  points. 

For  region  A,  roughness  lengths  were  calculated  separately  from  the 
vegetation  data  (SMR)  and  from  the  elevation  data  (TDR)  for  each  one  square 
kilometer  region  or  "plat".  About  12,000  square  kilometers  were  available 
with  both  elevation  and  vegetation  data.  In  region  B,  which  Includes  the 
smaller  region  A,  the  TDR  was  computed  for  each  square  kilometer.  Additional 
elevation  statistics  were  compiled  for  each  plat  In  both  regions  A  and  B, 
Including  the  average  elevation,  the  standard  deviation  of  the  set  of  eleva¬ 
tions  In  the  plat,  and  the  maximum  and  minimum  elevations  In  the  plat. 

The  grid  size  of  one  kilometer  Is  the  approximate  horizontal  scale  for 
boundary  layer  adjustment  to  changes  in  surface  roughness  to  occur.  Smaller 
"patches"  of  distributed  roughness  do  not  completely  alter  the  nature  of  the 
boundary  layer  as  does  a  larger  area  of  distributed  roughness  with  a  dimension 
of  one  kilometer  or  greater. 

All  results  for  both  regions  were  given  for  each  plat  as  Identified  In 
the  Universal  1 < «  sverse  Mercator  (UTM)  system  of  coordinates.  For  region  A 
this  format  presented  little  difficulty  since  the  Input  data  was  supplied  In 
the  same  coordinates.  Region  B  presented  more  of  a  problem.  This  area  was 
divided  Into  28  separate  sections.  The  boundaries  of  each  section  are 
constant  latitude  ahd  longitude  lines,  which  are  not  coincident  with  lines  of 
constant  UTM  coordinates.  However,  the  elevation  data  points  are  arranged 
along  UTM  coordinate  lines.  In  addition,  the  horizontal  distance  between 
adjoining  elevation  points,  63.5  m,  does  not  evenly  divide  into  1000m  (1  km). 
Because  of  this  arrangement,  some  square  kilometer  plats  lying  on  a  sectional 


boundary  may  not  have  a  full  set  of  elevation  heights.  In  some  cases  only  a 
few  data  points  were  available;  nevertheless,  the  TOR  and  elevation  statistics 
were  calculated  for  all  square  kilometers  with  at  least  three  good  data 
points.  The  number  of  plats  with  fewer  than  normal  number  of  data  points  was 
less  than  one  percent  of  the  total  number  of  plats. 

Because  of  the  horizontal  distance  between  elevation  points,  63.5  m, 
square  kilometers  are  Intersected  by  either  15  or  16  North-South  lines  of  data 
and  by  either  15  or  16  East-West  lines  of  data.  Thus  some  square  kilometers 
have  225  points,  some  240  points,  while  most  have  256  points.  The  number  of 
data  points  used  for  calculating  the  parameters  of  each  1  square  kilometer  was 
always  reported  for  each  plat. 

The  concept  of  a  terrain  derived  roughness  is  discussed  in  a  following 
section.  Although  a  variety  of  vegetation  based  roughness  measurements  have 
been  made,  there  Is  no  complete  set  of  roughness  lengths  that  corresponds  to 
many  commonly  encountered  types  of  vegetation.  It  Is  sometimes  possible  to 
categorize  a  specific  vegetation  form  or  type  for  which  no  roughness  measure- 


SURFACE  MORPHOLOGY 


The  smaller  area  data  base,  region  A,  has  31  codes  to  indicate  the  vege¬ 
tation  type  representative  of  each  grid  point.  The  roughness  lengths  asso¬ 
ciated  with  each  .egetation  code  were  obtained,  as  far  as  possible,  from 
published  data  and  assembled  in  a  table  of  roughness  lengths  by  vegetation 
code.  Published  values  of  the  roughness  length  of  vegetation  that  matched  a 
particular  vegetation  code  were  used;  many  of  the  values  are  listed  in  the 
ESDU  compilation  [5].  Several  vegetation  codes  had  no  corresponding  published 
value  for  roughness  length.  In  those  cases  the  appropriate  roughness  length 
was  estimated  by  combining  roughnesses  from  similar  vegetation  and  taking  an 
average. 

_2 

Agricultural  crops  were  listed  with  a  roughness  length  of  4.5x10  m. 

This  was  taken  to  represent  the  summer  value.  For  a  snow  covered  surface  in 

-3 

the  winter  the  value  was  2x10  m.  In  the  event  that  the  surface  was  not  snow 

_3 

covered  (based  on  climatology)  then  the  roughness  was  7.5x10  m.  Uncut 

grass,  z  *  2x10  m,  was  selected  as  representative  of  grassland,  pasture, 
o  -3 

and  meadows.  In  the  winter  with  snow  cover  the  z„  was  2x10  m  and  without 

-3  0 

snow  cover  It  was  7.5x10  m. 

Six  estimates  from  the  literature  of  roughness  length  for  coniferous 
forest  were  averaged  to  arrive  at  1.1  m.  Deciduous  forest  was  estimated  at 
1  m  for  the  summer.  The  ratio  of  summer  to  winter  roughnesses  of  "few  trees" 
(5.5  to  1)  was  used  to  estimate  the  fitter  roughness  length  of  a  deciduous 
forest  as  0.18  m.  Mixed  forest  was  assumed  to  be  50  percent  deciduous  and  50 
percent  coniferous  and  their  values  were  averaged  by  season.  Forest  clearings 
and  cutover  areas  were  assumed  to  be  50  percent  mixed  forest  and  50  percent 
Isolated  trees.  The  seasonal  values  are  simply  the  average  of  the  two. 

Orchards  were  considered  to  be  the  average  of  "few  trees"  and  "many 
-1  -2 

trees”  or  1.3x10  m  In  summer  and  1.4x10  m  In  the  winter.  The  value  for 
"many  hedges"  was  used  to  estimate  the  roughness  length  for  vineyards  and 
hop-gardens.  Brushland  and  scrubgrowth  (dense)  were  considered  to  be  the 


[5]  ESDU  7206,  1972,  "Characteristics  of  Wind  Speed  in  the  Lowest  Layers  of 
the  Atmosphere  near  the  Ground:  Strong  Winds,"  Eng.  Scl.  Data  Unit, 
Ltd.,  251  Regent  St.,  London,  WIR7AD. 


een  many  trees  and  hedges  and  many  hedges  or  1.4x10  1  m  in 

_2 

1/5.5)  of  this  value  in  winter  (2.5x10  m).  "Open  brushland"  was 

*  50  percent  brush  and  50  percent  grass.  The  resulting  seasonal 

i  average  of  the  two.  "Wetlands"  consisted  of  an  average  of  a  "few 

arge  expanses  of  water  or  Ice.  Values  for  the  description  “nearly 

_  3 

closely  spaced  low  growing  vegetation"  (9.5x10  m)  were  estimated 

between  "grassland/pasture/meadows"  and  "widely  spaced  low  growing 

"Abandoned  agriculture"  was  assumed  in  summer  to  be  long  grass 

-3 

ind  in  winter  to  be  cut  grass  (7.5x10  m).  "Peat  cuttings"  were 

>e  the  same  as  "grass/pasture"  (2x10  m).  Default  values  for  "no 

-2 

set  to  10  m.  The  roughness  table  summarizes  lengths  by  vegeta- 
>r  each  month  (Table  1).  In  transitional  seasons  (April-June  and 
ie  values  were  linearly  interpolated  between  respectively  March  and 
or  October  and  December  values.  This  procedure  was  used  for  all 
lants  and  trees. 

egion  A,  an  average  snow  cover  was  not  expected,  because  the  region 
ry  mean  tempeature  of  0°C  or  above  [6],  which  does  not  permit  snow 
/  long  on  the  ground. 

al  values  were  derived  from  the  vegetation  code  data.  We  calcu- 
nthly  characteristic  vegetation  roughnesses  for  each  square 
The  procedure  was  to  first  use  the  vegetation  roughness  lengths  to 
zQ  to  each  of  the  100  points  (region  A)  and  then  find  the  average 
An  ZQ' s,  called  L^.  The  "monthly  characteristic  vegetation 

zQ,  was  then  computed  as 

£o  *  exP(LAV)  * 

se  of  the  log  average  was  based  on  the  recommendation  by  Kung  [7] 
an  average  zQ  for  a  region  of  mixed  roughness  elements. 


ological  Office,  "Tables  of  Temperature,  Relative  Humidity, 

Itatlon  and  Sunshine  for  the  World,"  Part  III,  Europe  and  the 
,  London,  Her  Majesty's  Stationary  Office. 

E.  C.,  1963,  "Climatology  of  Aerodynamic  Roughness  Parameters  and 
Dissipation  In  the  Planetary  Boundary  Layer  of  the  Northern 
here.  Annual  Report,  Dept,  of  Meteorology,  University  of  Wisconsin, 
039-AMC -000878,  37-96. 
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TABLE  1 

ROUGHNESS  LENGTHS,  2 
( meters )  8 

MAR  APR 


APR 

MAY 

-3 

i 

7.5  x  10"3 

I  2.0  x  10' 

-3 

i 

7.5  x  10”3 

1.15  x  10' 

-2 

1.0  x  10"3 

|  2.5  x  10' 

0.18 

0.45 

0.64 

0.78 

.32 

0.40 

r2 

2.4  x  10"2 

5.8  x  10' 

,-2 

1.45  x  10"2 

3.6  x  10' 

,-2 

2.5  x  10"2 

6.3  x  10' 
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r2 
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1.25  x  10 

r4 
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4.5  x  10 

,-3 

1 
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,-3 
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Vegetation  Code 


TABLE  1  (Continued) 

ROUGHNESS  LENGTHS,  z 
(meters) 


4.S  x  10'2 

4.5  x  IQ"2 

4.5  x  10‘2 

4.5  x  10“2 

2.6  x  10“2 

7.5  x  10"3 

2.0  x  10'2 

2.0  x  10“2 

2.0  x  10-2 

2.0  x  10‘2 

1.4  x  10"2 

7.5  x  10“3 

>2 

5.5  x  10  c 

5.5  x  10-2 

5.5  x  10“2 

5.5  x  10“2 

3.25  x  10'2 

1.0  x  10"2 

1.1 

«» 

- 

_  -i 

- 

1.0 

1.0 

1.0 

1.0 

0.59 

0.18 

1.05 

1.05 

1.05 

1.05 

0.85 

0.64 

0.55 

0.55 

0.55 

0.55 

0.44 

0.32 

1.3  x  10"1 

1.3  x  10"1 

1.3  x  10*1 

1.3  x  10’1 

8.7  x  10’2 

2.4  x  10“2 

8.0  x  10"2 

8.0  x  10"2 

8.0  x  10‘2 

8.0  x  10‘2 

4.8  x  10“2 

1.45  x  10’2 

1.4  x  10"1 

1.4  x  10"1 

1.4  x  10'1 

1.4  x  10”1 

8.2  x  10‘2 

2.5  x  10*2 

8.0  x  10”2 

8.0  x  10”2 

8.0  x  10“2 

8.0  x  10’2 

4.8  x  10‘2 

1.6  x  10'2 

2.75  x  10"2 

2.75  x  10"2 

2.75  x  10“2 

2.75  x  10'2 

1.6  x  10~2 

0.5  x  10'2 

2.0  x  10"2 

- 

- 

- 

- 

-3 

1.0  x  10  J 

1.0  x  10"3 

1.0  x  10"3 

1.0  x  10"3 

5.9  x  10"4 

1.8  x  10"4 

9.5  x  10“3 

9.5  x  10’3 

9.5  x  10’3 

9.5  x  10’3 

5.6  x  10"3 

1.7  x  10"3 

5.00  x  10“2 

5.00  x  10"2 

5.00  x  10‘2 

5.00  x  10'2 

2.9  x  10"2 

7.5  x  10“3 

5.0  x  10~4 

- 

- 

- 

- 

4.0  x  10"1 

- 

- 

- 

- 

- 

5.0  x  10~4 

- 

- 

«k 

- 

m 

4.0  x  10"1 

m 

- 

m 

4 

5.5  x  10"1 

m 

I 

- 

m 

m 

8.5  x  10’1 

i 

I  «• 

- 

m 

- 

10*2 

- 

- 

w 

- 
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Seasonal  roughness  histograms  for  each  square  kilometer  were  also 
generated  for  seasons  defined  as:  Winter— Oecember,  January,  February; 
Spring— March,  April,  May;  Summer— June,  July,  August;  and  Fall— September, 
October,  November.  The  range  of  zQ  values,  which  comprise  six  orders  of 
magnitude,  were  divided  Into  12  bands.  For  each  season  the  total  number  of 
data  points  whose  vegetation  code  represented  a  roughness  length  within  each 
of  the  twelve  roughness  length  bands  was  determined.  A  histogram  was  then 
generated  of  roughness  length  band  versus  number  of  occurrences  at  data  points 
within  the  one  square  kilometer  for  a  three  month  period.  The  12  bands  for 
each  histogram  are  listed  In  Table  2,  and  Table  3  links  the  vegetation  code 
numbers  with  vegetation  description. 


As  an  example,  the  January  and  July  vegetation  roughness  calculations 
for  a  typical  square  kilometer  are  given  below.  The  vegetation  codes  for  this 


square  kilometer 

appear 

below  just  as  they  do 

on  the 

data 

tape. 

2 

2 

2 

2 

2 

2 

2 

2 

2 

11 

2 

2 

2 

2 

2 

2 

2 

11 

11 

11 

2 

2 

2 

2 

2 

2 

11 

11 

11 

11 

2 

2 

2 

11 

11 

11 

11 

11 

11 

11 

2 

11 

11 

11 

11 

11 

11 

11 

11 

11 

2 

11 

11 

11 

11 

11 

11 

11 

11 

11 

2 

2 

11 

11 

11 

11 

11 

11 

11 

11 

2 

2 

11 

11 

11 

11 

11 

11 

11 

11 

2 

11 

11 

11 

11 

11 

11 

11 

11 

11 

2 

2 

2 

2 

11 

11 

11 

11 

11 

11 

The  vegetation  consists 

of  grassland  (code  2) 

for  which  z 

0 

-  7.5x10" 

•3  , 

m  In 

2x10  in  in  summer  and  open  brushland  (code  11)  for  which 


winter  and  z„ 

-2° 

z  «  1.6x10  m  in  winter  and  z_ 

o  0  . 

dures  described  perilously, .  in  zQ  s  are  assigned  to  each  of  the  100  points. 

There  are  36  code  2  points  (in  z„  »  -4.89  winter,  in  z 

o  o 

code  11  points  (in  z^  *  -4.13  winter,  in  z 

0  0 

average  in  zQ  Is  derived  for  both  January  and  July  by  adding  the  100  separate 

in  z  and  dividing  by  100. 
o 


8.0x10  e  in  summer.  Following  the  proce- 

f  the  100  points. 
-3.91  summer)  and  64 
-2.53  summer).  Collectively  the 


The  "monthly  characteristic  vegetation  roughness",  zQ,  Is  then  deter 
mined  to  be  .01218  m  In  January  and  .04857  m  In  July. 
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TABLE  2 

VEGETATION  ROUGHNESS  BANOS 
Win  Max 


Band 

Min  z 

0 

Max  2 

0 

Un  z  ) 

0 

tn  z  ) 

0 

1 

10"5m 

5  x  10*5m 

-11.5 

-9.90 

2 

5  x  10“®m 

10‘4m 

-  9.90 

-9.21 

3 

10“4m 

.4 

5  x  10  m 

-  9.21 

-7.60 

4 

5  x  10”4m 

10"3m 

-  7.60 

-6.91 

5 

10"  3m 

5  x  10  3ffl 

-  6.91 

-5.30 

6 

5  x  I0"3m 

,A-2 

10  m 

-  5.30 

-4.61 

7 

10"^m 

-2 

5  x  10  m 

-  4.61 

-3.00 

8 

5  x  10"2ra 

10’1® 

-  3.00 

-2.30 

9 

lO'1® 

5  x  10 

-  2.30 

-0.693 

10 

5  x  10" ^m 

1  m 

-  0.693 

0 

11 

1  m 

5  m 

0 

1.61 

12 

5  m 

10  m 

1.61 

2.30 

Typical 

Features 

Ice,  Mud  Flats 

Calm  Sea 

Snow  Surface 

Grass  Crops 

Crops 

Forest 
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TABLE  3 

VEGETATION  DESCKIPTION  BY  CODE 


Description  Code 


No  data 

0 

Agriculture,  cropland 

1 

Grassland,  pasture,  meadows 

2 

Grassland,  scattered  trees 

3 

Coni  fero’js  forest 

4 

Deciduous  forest 

5 

Nixed  forest 

6 

Forest  clearings,  cutover  areas 

7 

Orchards 

8 

Vineyards,  hop-gardens 

9 

Brushland,  scrub  growth  (dense) 

10 

Brushland,  scrub  growth  (open) 

11 

Wetl ands 

12 

Peat  cuttings 

13 

Nearly  barren  w/wldely  spaced 

low  growing  vegetation 

14 

Nearly  barren  w/closely  spaced 

low  growing  vegetation 

15 

Abandoned  agriculture 

16 

Bare  ground,  sand  dunes 

17 

/ 

18 

19 

20 

21 

<  Undefined 

22 

23 

24 

25 

< 

26 

Surface  mines 

27 

Open  water 

28 

Villages 

29 

Towns 

30 

Cities 

31 
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TERRAIN  DERIVED  ROUGHNESS 

The  computation  of  roughness  lengths  from  elevation  data  1$  based  upon  a 
formula  given  by  Lettau  [4]  for  estimating  the  roughness  for  distributions  of 
more-or-less  Identical  elements.  The  relationship  Is 

(1)  *  0.5h  S/A 

o 

where 

h  Is  the  average  vertical  extent  (height)  of  an  element  (m), 

S  Is  the  average  projected  area  of  an  element  perpendicular  to  the  wind 
2 

(m  ),  and 

2 

A  Is  the  average  horizontal  ("lot")  area  (m  )  occupied  by  an  element. 

As  Lettau  points  out.  equation  (1)  Incorporates  not  only  the  height  of  the 
elements  but  also  a  measure  of  their  presented  area  and  their  spacing  density. 
This  formula  should  then  be  preferred  over  the  exponential  relationships  (e.g. 
zq  -  ah**)  given  by  Kung  [7]  and  others  which  Involve  only  the  element  height. 

In  our  adaptation  of  Lettau* s  formula  the  elevation  data  Is  first  pro¬ 
cessed  to  obtain  estimates  of  h,  S,  and  A  over  a  specified  area.  The  data  Is 
supplied  as  elevations  In  meters  above  sea  level  given  on  a  square  grid.  The 
spacing  between  grid  points  Is  i  meters  and  the  lines  formed  by  the  rows  of 
points  run  generally  North-South  and  East-West.  Calculations  to  compute  a 

z  are  performed  on  each  row  (South  to  North  traverse)  and  then  on  each  column 
o 

(West  to  East  traverse)  of  the  elevation  data.  The  20  to  32  values  (the  exact 

number  depends  on  the  grid  Interval)  are  averaged  to  produce  the  "typical" 

z  for  the  square  kilometer.  This  cross  pattern  is  adopted  In  order  to  pro- 
o 

duce  an  average  estimate  that  Is  Independent,  as  much  as  possible,  of  wind 
direction.  The  expressions  below  which  apply  to  a  single  traverse  are  not 
dependent  upon  which  direction,  e.g.  South  to  North  or  North  to  South,  the 
traverse  Is  made. 

To  estimate  h  we  first  compute  the  number  of  "peaks  and  valleys",  N, 
encountered  In  a  traverse.  The  first  elevation  In  a  traverse  Is  considered  a 
peak  If  the  second  elevation  Is  lower  and  a  valley  If  the  second  Is  higher. 

For  the  second,  third,  fourth,  etc.  points  In  a  traverse,  a  peak  Is  counted 
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when  both  the  preceding  and  following  points  around  a  given  point  are  lower  In 
elevation.  A  valley  Is  counted  when  both  neighboring  points  are  higher.  In 
any  other  situation  no  increase  Is  made  In  N.  The  final  point  In  a  traverse 
Is  treated  In  the  same  way  as  the  first  point.  It  Is  a  peak  If  It  Is  higher 
than  the  next  to  last  point  and  a  valley  If  It  Is  lower.  As  the  value  of  N  Is 
computed  for  each  traverse  a  running  sum  of  the  total  vertical  excursion,  E, 

Is  kept.  E  Is  the  total  upward  and  downward  distance  In  meters  covered  In 
moving  along  the  row  of  elevation  points.  The  estimate  of  h  Is  then  computed 
as 

(2)  "-anblT 

-2 

When  N  Is  one  or  zero  we  have  a  flat  traverse  and  a  default  value  of  1x10  m 
Is  set  for  zQ.  the  remaining  computations  for  the  traverse  being  Ignored. 

For  the  average  projected  area  of  the  terrain  "bumps"  we  multiply  E/2  by 
the  width  of  the  grid  spacing,  t,  and  divide  by  N-l. 

I3»  S  ■  I  •  (H-l) 

We  are  thus  assuming  the  bumps  to  be  two-dimensional  "washboard".  The  lot 
area  per  bump  Is  the  traverse  horizontal  area  divided  by  (N-l),  or 

(4)  A  »  U* lkm) /(N-l) 

Combining  the  three  numbers  In  Lettau's  formula,  eq  (1),  we  obtain  the 
TOR  for  the  traverse  In  meters 

(5)  zQ  *  E2/[8000(N-1)]  . 

For  a  square  kilometer  the  reported  TOR  Is  the  average  zQ  for  the  M  (20  to 
32)  traverses, 
or 

<*'  s  •  ijj'Vi  • 

As  a  further  Indication  of  the  roughness  of  the  terrain  we  compute  the 
average  number  of  peaks  and  valleys,  N,  from  the  values  for  each  traverse 


-1 M 

N  a  «r  £  (N)j  • 

"l«l  1 

This  value  was  Included  with  each  calculation  of  the  TDR. 

As  an  example  of  the  above  procedure  consider  the  following  100  eleva 
tlons  from  a  typical  one  square  kilometer  plat: 


160 

153 

147 

146 

151 

141 

139 

140 

139 

150 

151 

148 

142 

140 

145 

139 

136 

138 

143 

150 

140 

140 

136 

134 

135 

135 

139 

142 

143 

149 

134 

134 

133 

134 

135 

141 

142 

149 

156 

165 

133 

135 

135 

138 

150 

161 

152 

167 

180 

183 

138 

138 

140 

150 

177 

180 

170 

181 

184 

184 

149 

143 

148 

154 

180 

181 

184 

184 

184 

180 

155 

150 

155 

159 

171 

178 

183 

184 

181 

169 

160 

162 

163 

168 

168 

175 

182 

184 

180 

162 

160 

166 

171 

175 

178 

180 

185 

182 

169 

158 

is  case 

there 

i  are  2C 

l  traverses,  10 

South- 

North 

and  10 

East-West. 

number  the  traverses  from  the  left  for  South-North  and  from  the  bottom  for 
East-West,  the  values  of  N,  E,  and  zQ  computed  by  equations  (2)  through  (5) 


Traverse 


E  (meters) 


r  (meters) 


i. 


(E-W) 


1 

3 

52 

.1690 

2 

3 

46 

.1323 

3 

4 

61 

.1550 

- 

4 

3 

51 

.1626 

5 

2 

66 

.5445 

% 

6 

4 

68 

.1927 

7 

2 

33 

.1361 

8 

2 

21 

.0551 

9 

5 

39 

.0475 

10 

7 

42 

.0268 

Finally  the  characteristic  roughness  length  for  this  plat  is  the  mean  of  these 
20  values 

z  *  .2015  meters, 
o 

and  the  average  number  of  peaks  and  valleys  is 
N  *  3.35  . 

The  additional  statistics  for  this  plat  are:  mean  elevation  Is 
156.9  meters,  the  standard  deviation  Is  17.53  meters,  and  the  maximum  and 
minimum  elevations  are  185  and  133  meters  respectively. 
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DISCUSSION 


Calculations  of  Surface  Morphology  Roughnesses  (SMR)  are  consistent  with 
expectation.  Without  regard  to  season  or  location,  SMR  lengths  range  predomi¬ 
nately  between  .01  and  1  meters  with  few  exceptions.  They  are  greatest  in 
summer  months  for  cities  and  forest  areas  and  least  in  winter  months  for  grass 
covered  plats. 

Terrain  derived  roughnesses  (TOR)  are  also  consistent  with  expectations. 
They  are  greatest  for  hilly  country  and  least  for  flat  areas.  Their  range  is 
from  .01  meter  to  above  1  meter.  It  is  interesting  to  compare  the  TDR's  to 
the  more  traditional  SMR's  especially  In  hilly  country.  The  TDR's  are  con¬ 
siderably  larger.  It  Is  reasonable  to  expect  that  the  vegetative  roughness  in 
hilly  country  Is  of  secondary  Importance. 

TDR's  for  region  A  with  a  granularity  of  100  meters  are  generally  larger 
than  those  calculated  for  Identical  square  kilometer  areas  using  region  B  data 
whose  granularity  Is  63.5  m.  In  addition,  the  standard  deviations  of  eleva¬ 
tion  height  using  100  meter  data  Is  also  larger  than  that  calculated  using 
63.5  meter  data.  This  comparison  suggests  that  both  the  TDR  and  elevation 
standard  deviation  are  sensitive  to  the  Interval  between  data  points.  Further 
study  of  this  sensitivity  to  granularity  Is  required. 
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APPENDIX 

COMPUTER  PROGRAM  LISTINGS 

This  appendix  contains  the  FORTRAN  listings  of  the  computer  programs 
used  to  compute  the  roughness  lengths  and  other  terrain  data.  Program  VEGRUF 
was  used  for  region  A  which  contained  both  elevation  and  surface  type  data. 
Program  TDRUFF  was  used  for  region  B  which  contained  only  elevations.  Both 
programs  were  developed  and  executed  on  a  Digital  Equipment  Company  VAX/11-780 
computer  using  the  VAX/VMS  Version  2  operating  system  and  VAX  FORTRAN-IV  PLUS. 
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100 
200 
300 
310 
320 
330 
340 
390 
360 
370 
380 
390 
400 
500 
600 
700 
800 
810 
820 
830 
900 
1000 
1  100 
1200 
1300 
1400 
1500 
1600 
1700 
1800 
2000 
3500 
3600 
3650 
3700 

3709 
3706 

3710 
3720 
3730 
3740 
3780 
3800 
3900 
4000 
4100 
4200 
4300 
4400 
4500 
4600 
4700 
4800 
4900 
9000 
5100 
9200 
9230 
9300 
5400- 
5500 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


PROGRAM  VECRUF 

VECRUF  COMPUTES  VALUES  OP  SURFACE  MORPHOLOGY 
ROUGHNESS.  TERRAIN  DERIVED  ROUGHNESS.  AND  TERRAIN 
STATISTICS  FOR  ELEVATION  AND  SURFACE  CHARACTERISTICS 
DATA  SUPPLIED  ON  A  SOUARE  GRID 

VECRUF  MAS  CREATED  FOR  THE  US  ARMY  ATMOSPHERIC 
SCIENCES  LABORATORY  BY  THE  UNIVERSITY  OF  DAYTON 
RESEARCH  INSTITUTE  31  MARCH  1981 

NOTE:  SOME  PARTS  OF  THIS  CODE  MAY  BE  SPECIFIC  TO 
DEC  VAX/VMS  FORTRAN. 


INPUT: 

UNIT  1:  RAW  DATA  TAPE  WITH  DATA  SPACED  AT  100  METER 
GRANULARITY 

DATA  IS  TO  BE  SUPPLIED  AS  97  KM  LONG  SCAN  LINES 
ONE  LINE  PER  RECORD 
UNIT  9:  VSER  IDENTIFICATION  OF  RUN 
RECORD  1: 

IUTMX-UTM  X  COORDINATE  OF  THE  SOUTHWEST  CORNER  OF  THE  DATA  AREA 
IUTMY-UTM  Y  COORDINATE  OF  THE  SOUTHWEST  CORNER  OF  THE  DATA  AREA 

OUTPUT: 

UNIT  2:  ROUGHNESS  DATA-  OUTPUT  FOR  EACH  SOUARE  KM. 

UNIT  6:  ECHO  OF  USER  IDENTIFICATION  ABOVE 


COMMON/COLUMN/CR 1D<970» 10. 2). NPTINC.  NR C INC.  NPTKM.  NRCKH.  XNTOT 
COMMON /KOUT/HMEAN( 97 > » HSIGMA(97>.  HMAX (97).  HMIN(97>»  HZER0C97), 

1  HNPV (97) 

COMMON/RUFOUT/IVC<97).NU<97>.  ZHAT <12.  97),  IHISTO< 12.  4.  97) 
DIMENSION  MEAN < 97 ) , ISIGMAC97),  IHMAX(97>,  IHM1N<97),  IHZER0<97>. 

1  I ZHAT < 12,  97),  I HNPV (97) 

INTEGER *2  GRID 

EQUIVALENCE  ( HME AN, MEAN ) .  (HSICMA.  IS1GMA),  (HMAX,  I HMAX ) , 

1  (HMIN, IHMIN).  (HZERO,  IHZERO).  (ZHAT. IZHAT). 

2  (HNPV, I HNPV) 

IFlLER-lIllllllll 

NPTINC-1 

NRCINC-1 
NPTKM- 10/NPT1NC 
NRCKM-10/NRCINC  , 

XNTOT-FLOAT (NPTKM*NRCKM) 

0PEN(UNIT-2.  BLOCKSI ZE-20400.  RECORDSI ZE-20400.  RECORDTYPE- 'FIXED 
1  CARR I ACECONTROL- 'NONE '. STATUS- 'NEW' ) 

READ  IN  UTM  COORDINATES  OF  THE  SOUTHWEST  CORNER  AND 
ECHO  ON  OUTPUT. 

READ( 5,  10) IUTMX.  1UTMY 
10  FORMAT (2110) 

WRITE(6. 20) IUTMX.  IUTMY 

20  FORMAT ( 1H1,  'SW  UTM  COORDINATES  -  (M10.  M10.  ')') 

WRITE*  2, 1100) (IUTMX.  IUTMY.  1-1.  1020) 

READ  IN  AMD  COMPUTE  THE  MEAN  HEIGHT.  AND  ROUGHNESS  HEIGHT  FOR 
EACH  SOUARE  KM  OF  125  97KM.  COLUMNS. 
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•*  t*  n 


DO  1000  T-l.  126 


SUBR  GETCOL  ARRANGES  THE  SCAN  LINES  INTO  A  1  KM  WIDE  !E-W>  COLUMN 
OF  DATA  97  KM  HIGH  (N-S).  DATA  IS  STORED  IN  ARRAY  'GRID' 

CALL  GETCOL 

IF < 1  LT.  120)00  TO  1000 

SUBR  HCOMP  COMPUTES  THE  OUTPUT  DATA  OBTAINED  FROM  THE 
ELEVATION  DATA.  I.  E.  THE  TDR  AND  TERRAIN  STATISTICS. 

CALL  HCOMO 

WRITE  <6, 40)HMEAN.  HSICMA.  HMAX.  HMIN.  HZERO.  HNPV 
FORMAT ( 1 X.  1 OF 1 1 .  4) 

SUBR  CETRUF  COMPUTES  THE  OUTPUT  DATA  OBTAINED  FROM 
THE  SURFACE  TYPE  CODES.  I.E.  THE  SMR 

CALL  CETRUF 

WRITE! 6. SO) IVC.  NU.  ZHAT.  IHISTO 

FORMAT! 19 ( IX.  lOI 10/ ) .  1 1 X, 4 1 10/ ) , 97 < 12F 1 1 .  4/ ) .  < 1 X .  1 21 1 0/ ) ) 

SCALE  THE  OUTPUT  VARIABLE  AND  CONVERT  TO  INTEGER 
THIS  IS  A  SPACE  SAVING  PROCEDURE  FOR  THE  OUTPUT  TAPE 

DO  150  J-1.97 

MEAN! J)-IFIX!HMEAN!U)*10.  ♦.  5) 

1SICMA! J)=IFIX!HSICMA! J)*100.  +  5) 

IHMAX!U)-1FIX!HMAX! J)+.  1) 

I  HMIN!  J)  -a  IFI X  !  HMIN!  J)  ♦.  1) 

IHNPV! J)-IFIX!HNPV!J)*100.  ♦.  5) 

IHZERO ! J ) “IF I X ! HZERO! J)*l.  E6«v  5) 

DO  140  K-l.  12 

1ZHAT!K.  J)“IFIXIZHAT!K.  J>*1.  E6«-.  5) 

\0  CONTINUE 
>0  CONTINUE 
NSTART  *—29 
NEND-0 
DO  160  M»l,3 
NSTART "NSTART *30 
NEND-NEND+30 

WRITE <2.  1100)! <MEAN!J),  ISICMA(J),  IHMAX!J),  IHMIN! J),  IHNPV!J). 

1  IHZERO! J),  IVC! J). NU! J).  ! I ZHAT !K, J).  K"l.  12). 

2  ! ! IHIST01K,  L.  J).  K“1 .  1 2).  L»l.  4) ) .  U-NSTART,  NEND) 

>0  CONTINUE 

WRITE! 2.  1100) ! I MEAN! J).  ISIGMA!U>,  IHMAX!U).  IHMIN! J).  IHNPV! J). 

IHZERO!U).  IVC! J). NU! J).  ! I ZHAT !K. J  . K-l,  12). 

! ! IHIST01K,  L.  J).K»1.  12).  L-l.  4) ), J-91. 97). 

! IFILER. N-l, 1564) 

K)  CONTINUE 
•O  FORMAT !2C40I 10) 

STOP 

END 
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100 
200 
300 
400 
300 
600 
700 
800 
900 
1000 
1 100 
1200 
1300 
1400 
1500 
1600 
1700 
1800 
1900 
2000 
2100 
2200 
2300 
2400 
2500 
2600 
2700 
2800 
2900 
3000 
3100 
3200 
3300 
3400 
3300 
3600 
3610 
3620 
3630 
3640 
3700 
3800 
3900 
4000 
4100 
4200 
4210 
4500 
4600 
4700 
4800 
4900 
3000 
3100 
3200 
3300 
3400 
3300 
3600 
3700 
3800 


SUBROUT  1NE  CETCOL 
C 

C  THIS  SUROUTINE  FILLS  THE  ARRAY  CRIB  WITH  THE  HEIGHT  AND 

$  ROUGHNESS  TYPE  FOR  EACH  DESIRED  DATA  PT.  IN  A  COLUMN  CF  97  SQUARE  KMS 

C 

C  INPUT  AND  OUTPUT  INCLUDED  IN  COMMON/COLUMN/ 

C 

C  INPUT: 

C  NPTINC-THE  INCREMENT  ADDED  TO  THE  POINT  COUNTER  TO  DETERMINE 

C  WHICH  POINTS  WITHIN  A  RECORD  WILL  BE  USED. 

C  NRCINC-THE  INCREMENT  ADDED  TO  THE  RECORD  COUNTER  TO  DETERMINE 

C  WHICH  RECORDS  WILL  BE  USED. 

C  NPTKM  -THE  RESULTANT  NUMBER  OF  POINTS  TO  BE  EXTRACTED  PER  KM 

C  FROM  A  RECORD. 

C 

C  OUTPUT: 

C  CRIDO- 

C  THE  FIRST  SUBSCRIPT  VARIES  FROM  1  TO  NPTKM*97  AND  INDEXES 

C  THE  POINTS  OF  A  97  KM  SPAN  IN  THE  N-S  DIRECTION. 

C  THE  SECOND  SUBSCRIPT  VARIES  FROM  1  TO  10/NSClNC  AND  INDEXES 

C  THE  RECORDS  WITHIN  A  KM  SPAN  IN  THE  E-W  DIRECTION 

C  THE  THIRD  SUBSCRIPT  DENOTES  THE  TYPE  OF  DATA  STORED  IN  GRID. 

C  1  IMPLIES  ELEVATION  IN  METERS  ABOVE  MEAN  SEA  LEVEL. 

C  2  IMPLIES  SURFACE  TYPE. 

C 

COMMON /COLUMN/ CR 1 D  < 970.  10,  2),  NPTINC,  NRCINC.  NPTKM,  NRCKM, XNTOT 
INTEGER *2  IBUFU942) 

BYTE  BUFFER <3884) 

BYTE  IEL6YT<2>.  ISURBY<2> 

EQUIVALENCE ( 1BUF ( 1 ) , BUFFER < 1 ) ) ,  < IELBYT 1 1 ) ,  IELEV),  (ISURBYCl),  ISURF) 

INTEGER *2  IELEV. ISURF, GRID 

PARAMETER  I0«  READLBLK  -  '21 'X 

PARAMETER  EOF- '870 'X 

PARAMETER  BUFS1 ZE-3884 

INTEGER *2  CHANNEL. I0SB<4) 

INTEGER *4  SYS9ASS1GN,  SYS9QI0U 
LOGICAL  L9ASSI0N 

DATA  L9ASSIGN/.  TRUE.  /. NREC3/0/ ,  NXTREC/ 1 / 

I F  <  L« ASS I GN ) THEN 

L9ASSIGN-.  FALSE. 

IRET«SYS*ASSIGN< 'F0R001 '.CHANNEL.  .  > 

IF  (IRET.  NE.  1.  THEN 
WRITE (6.  3 > I RET 

5  FORMAT < '  SYS«ASSIQN  ERROR'.  ZB) 

ST  37 
ENOIF 
ENDIF 

J-0 

DO  100  I-l. 10 
C 

C  READ  INPUT  RECORD.  < THERE  ARE  10  PER  KM> 

C 

10  I  RET  «SYS*Q  1 OW  < ,  7.VAL  <  CHANNEL  >.  7.VAL  <  10*_READLBLK>.  IOSB.  .  . 

1  XREF < BUFFER >.  XVAL < BUFSI ZE >....  > 

IF< IRET.  WE.  1)THEN 
WRITE (6. 13) IRET 

13  FORMAT! '  READ  OlOW  ERROR'.  ZB) 

STOP 

ENDIF 

IFtlOSBU).  EQ  EOF.  OR.  lOSB(S).  EQ  0 ) THEN 
CLOSE < UNIT- 1 ) 
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WRITE <6. *>NRECS.  RECORDS  READ 

STOP 

END  IF 

NRECS-NRECS*1 

IFCNXTREC.  NE.  NRECSiGO  TO  100 

NXTREC-NXTREC+NRCINC 

IF < NRECS  LE  1190)00  TO  100 

J-J+l 

INDEX-0 

LOOP  OVER  97  KHS  WITHIN  RECORD 
DO  50  K-1.97 

<10  PTS  PER  KM*4  BYTES  PER  PT)-40  BYTES  PER  KM 
DO  40  L-l. 40.  NPTINC*4 
NXTPNT- <K-1>*40*L 
I ELBYT  <  2 ) -BUFFER  < NXTPNT ) 

IELBYT < 1 '-BUFFER <NXTPNT+1 ) 

I SURB Y  <  2 ' -BUFFER  <  NXTPNT *2 ) 

I SURBY  < 1 > -BUFFER  < NXTPNT+3 > 

INDEX- INDEXfl 

OR ID< INDEX.  J.  11-IELEV 

CRID< INDEX.  J.2>-1SHFTUSURF, -11) 

40  CONTINUE 
50  CONTINUE 
100  CONTINUE 
RETURN 


100 

200  c 

300  C 
400  C 
500  C 
600  C 
700  C 
800  C 
900  C 
1000  C 
1100  C 
1200  C 
1300  C 
1400  C 
1500  C 
1600  C 
1700  C 
1800  C 
1900  C 
2000  C 
2100  C 
2200  C 
2300  C 
2400  C 
2500  C 
2600  C 
2700  C 
2800  C 
2900  C 
3000  C 
3100  C 
3150  C 
3200  C 
3300 
3400 
3450 
3500 
3600 
3700 
3800 
3900 
4000 
4050 
4100 
4200 
4300 
4400 
4500 
4600 
4700 
4800 
4900 
5000 
5100 
5200 
5300 
5400 
5500 
5600 
5700 
5800 


SUBROUTINE  HCOMP 

THIS  SUBROUTINE  COMPUTES.  THE  OUTPUT  VARIABLES  DERIVED  FROM  TH 
RAW  ELEVATION  DATA  FOR  EACH  OF  THE  97  SQUARE  KMS  CONTAINED  IN 
ARRAY  GRID 


INPUT: 

NPTKM-THE  NUMBER  OF  POINTS  IN  THE  N-S  DIRECTION  PER  HM 
NRCKM-THE  NUMBER  OF  POINTS  IN  THE  E-W  DIRECTION  PER  KM 
XNTOT-THE  TOTAL  NUMBER  OF  POINTS  PER  KM 
GRID!  >-" 

THE  FIRST  SUBSCRIPT  VARIES  FROM  1  TO  NPTKM*97  AND  INDEXES 
THE  POINTS  OF  A  97  KM  SPAN  IN  THE  N-S  DIRECTION. 

THE  SECOND  SUBSCRIPT  VARIES  FROM  1  TO  10/NRCINC  AND  INDEXES 
THE  RECORDS  WITHIN  A  KM  SPAN  IN  THE  E-W  DIRECTION.  . 

THE  THIRD  SUBSCRIPT  DENOTES  THE  TYPE  OF  DATA  STORED  IN  GRID 

1  IMPLIES  ELEVATION  IN  METERS  ABOVE  MEAN  SEA  LEVEL. 

2  IMPLIES  SURFACE  TYPE.  (NOT  USED  IN  THIS  SUBROUTINE) 


OUTPUT: 

N.  B.  TH?  LAST  SUBSCRIPT  OF  EACH  OF  THE  OUTPUT  ARRAYS  VARIES 
FROM  1  TO  97  WITH  ONE  VALUE  FOR  EACH  OF  THE  97  SQUARE  KMS 
FROM  SOUTH  TO  NORTH  RESPECTIVELY. 

HMEAN ( 97 ) “THE  MEAN  ELEVATION  (METERS) 

HSICHA(97)»THE  ELEVATION  STANDARD  DEVIATION  (METERS) 

HMAX ( 97 ) “THE  MAXIMUM  ELEVATION  (METERS) 

HMIN( 97 ) “THE  MINIMUM  ELEVATION  (METERS) 

HZER0(?7)»THE  TERRAIN  ROUGHNESS  PARAMETER  (METERS) 
HNPV(97)-THE  AVERAGE  NUMBER  OF  PEAKS  AND  VALLEYS/SQ.KM 

COMMON/COLUMN/ OR ID<970. 10.  2).  NPTINC.  NRCINC.  NPTKM,  NPCKM.  XNTOT 
COMMON/HOUT/HMEAN ( 97 ) . HSI CMA ( 97  > .  HMAX ( 97 ) . HMI N ( 97 ) .  HZERO ( 97 ) . 

HNPV ( 97) 

INTEGERS  GRID 
DIMENSION  HTRAV(20> 

DO  lOOO  U-1.97 
HMEAN< J)«0.  0 
HSICMA(J)-0.  0 
HZERO( J)-0.  0 
NPVDIR-0 

HMAX (J> -FLOAT (OR ID ( (U-l ) *NPTKM+I.  1.  I)) 

HMIN( J)-HMAX(J) 

DO  90  L-l.NRCKM 
DT0TAL*0.  0 
NPNV-0 

HNEXT-FL0AT(GRID((U-1)*NPTKM*I,L.  1>> 

DO  80  •*.*!.  NPTKM 
H-HNEXT 

HMEAN  ( J  )  -HMEAN  (  U  >  «+( 

HSICMA( J)«HSIGMA( J)+H**2 
IF (H.  GT.  HMAX(U) )HMAX( J)-H 
1F(H.  LT.  HMIN(J))  HMIN(J)-H 
IF(K.  EQ.  1)  THEN 

HNEXT-FL0AT(CR1D< (U-i ) *NPTHM*H*1>  L.  1>> 

DIFF*HNEXT-H 
IF(DIFF.  NE.  0.  0)THEN 
NPNV-NPNV*! 

DTOT AL-DTOT AL ♦ ABS ( D I FF ) 
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5900 

END  IP 

6000 

HrREV-H 

6  tOO 

ELSE  IriK.  LT  NPTKH)  THEN 

6200 

HweXT-FLOAT  (GRID!  1  J-l  )*NPTKH«-K*1>  L«  1>> 

6300 

1K<  <H.  CT  HPREV.  AND  H  CT.  HNEXT).  OR. 

6400 

1  <H,  LT.  HPREV  AND.  H.  LT.  HNEXT)  )  NPNV-NPNV* 1 

6500 

DIFF-HNEXT-H 

6600 

DtQTAL-DT0TAL*A8S<  DIFF  > 

6700 

HPREV-H 

6900 

ELSE 

6900 

DIFF-H-HPREV 

7000 

IF < D1FF.  NE.  0.  OJNPNV-NPNV+l 

7100 

END  IF 

7200 

80 

CONTINUE 

7300 

IF  (NPNV.  EQ.  0.  OR.  NPNV.  EQ.  1  >  THEN 

7400 

HTRAV<L>-1.  OE-2 

7500 

ELSE 

7600 

HTRAV<L)-DT0TAL**2/< CNPNV-1 >*8000.  > 

7700 

END  IF 

7750 

NPVD 1R-NPVD I R+NPNV 

7  BOO 

1F<HTRAV(L).  EQ.  0.  0>HTRAV(L>-1.  OE-2 

7900 

90 

CONTINUE 

8000 

HS 1 CMA <  J ) “SORT < ( HSIGMAIJ)  —  (HHEAN<  J) **2/XNTQT ) ) / t  XNTOT— 1 .  0) ) 

8100 

HUE AN ( J ) -HUE AN <  J ) / XNTOT 

8200 

DO  190  K-l.  NPTKH 

8300 

DTOTAL-0.  0 

8400 

NPNV-0 

8500 

HNEXT-FL0AT<CR1D< (J-l )*NPTKH+H.  1.1>) 

8600 

DO  180  L-l.NRCHH 

8700 

H=*  HNEXT 

3800 

IF<L.  EQ.  1>  THEN 

8900 

HNEXT-FLOAT (GRID! < J-l ) *NPTHH+K.  L«- 1 .  1  >  ) 

9000 

DIFF-HNEXT-H 

9100 

IF1DIFF.  NE.  0.  0>THEN 

9200 

NPNV-NPNV* 1 

9300 

DTOTAL-DTOTAL+ABS(DIFF> 

9400 

END  IF 

9500 

HPREV-H 

9600 

ELSE  I  =  :L.  LT.  NRCHH)  THEN 

9700 

HNEXT-FLOAT  (GRID!  ( J-l  >*NPTKH*H.  L«-l.  1)> 

9800 

IF (  (H.  GT.  HPREV.  AND.  H!  CT.  HNEXT).  OR. 

9900 

1  (H.  LT.  HPREV  AND.  H.  LT.  HNEXT)  )  NPNV-NPNV-M 

lOOOO 

DIFF-HNEXT-H 

10100 

DTOT AL-DTOT AL+ABS <  D I FF ) 

10200 

HPREV-H 

10300 

ELSE 

10400 

DIFF-H-HPREV 

10500 

IF1DIFF.  NE.  0.  0)  NPNV-NPNVM 

10600 

END  IF 

10700 

180 

CONTINUE 

10800 

IF < NPNV.  EQ.  0.  OR.  NPNV.  EQ.  1)  THEN 

10900 

HTRAV<NRCKH«-K>-1.  OE-2 

11000 

ELSE 

11100 

HTRAV<NRCKH*K)-DT0TAL«*2/» < NPNV- 1 >*8000.  ) 

11200 

END  IF 

11250 

NPVDIR“NPVDIR*NPNV 

11300 

IF<HTRAV<NRCKH*H).  EQ.  0.  0>HTRAV<NRCKH*K>-1.  OE-2 

11400 

190 

CONTINUE 

11500 

DO  300  K-l.  NRCKH+NPTHH 

11600 

HZEROC J)-HZERO< J) +HTRAV1H) 

11700 

300 

CONTINUE 
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1 1600 
11850 
11900 
12000 
12100 


HZERO  <  J ) -HZERO <  J > / < NRCKrt»NRTKH > 

HNPV  <  Jl-FLOAT  <  NPVDIR  >  /FLOAT  ( NRCKHfNPTKHl 
1000  CONTINUE 
RETURN 
END 
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100 

200  C 
300  C 
400  C 
500  C 
600  C 
700  C 
800  C 
900  C 
1000  C 
1 100  C 
1200  C 
1300  C 
1400  C 
1500  C 
1600  C 
1700  C 
1800  C 
1900  C 
2000  C 
2100  C 
2200  C 
2300  C 
2400  C 
2500  C 
2600  C 
2700  C 
2800  C 
2900  C 
3000  C 
3100  C 
3200  C 
3300  C 
3400  C 
3500  C 
3600  C 
3700  C 
3800  C 
3900  C 
4000  C 
4100  C 
4200- 
4300 
4400 
4500 
4600 
4700 
4800 
4900 
5000 
5100 
5200 
5300 
5400 
5450 
5500 
5600 
5700 
5800 
5900. 
5950 


SUBROUTINE  CETRUF 

THIS  SUBROUTINE  COMPUTES  THE  OUTPUT  VARIABLES  DERIVED  FROM  THE 
SURFACE  TYPE  DATA  FOR  EACH  OF  THE  97  SQUARE  KMS  CONTAINED 
IN  ARRAY  GRID. 


INPUT: 

NPTKM-THE  NUMBER  OF  POINTS  IN  THE  N-6  DIRECTION  PER  KM 
NRCKM-THE  NUMBER  OF  POINTS  IN  THE  E-W  DIRECTION  PER  KM 
6R1D<  >- 

THE  FIRST  SUBSCRIPT  VARIES  FROM  1  TO  NPTKM«97  AND  INDEXES 
THE  POINTS  OF  A  97  KM  SPAN  IN  THE  N-S  DIRECTION. 

THE  SECOND  SUBSCRIPT  VARIES  FROM  1  TO  10/NRCINC  AND  INDEXES 
THE  RECORDS  WITHIN  A  KM  SPAN  IN  THE  E-U  DIRECTION 
THE  THIRD  SUBSCRIPT  DENOTES  THE  TYPE  OF  DATA  STORED  IN  GRID 

1  IMPLIES  ELEVATION  IN  METERS  ABOVE  MEAN  SEA  LEVEL. 

(NOT  USED  IN  THIS  SUBROUTINE) 

2  IMPLIES  SURFACE  TYPE. 

OUTPUT: 

N.  B.  THE  LAST  SUBSCRIPT  OF  EACH  OF  THE  OUTPUT  ARRAYS  VARIES 
FROM  1  TO  97  WITH  ONE  VALUE  FOR  EACH  OF  THE  97  SQUARE  KMS 
FROM  SOUTH  TO  NORTH  RESPECTIVELY. 

I VC (97) -THE  CHARACTERISTIC  SURFACE  TYPE.  (MOST  COMMON) 
NU(97)-THE  NUMBER  OF  UNDEFINED  SURVACE  TYPE  CODES. 

ZHAT(12. 97)-THE  SET  OF  MONTHLY  CHARACTERISTIC  SURFACE 
ROUGHNESS  VALUES.  (ONE  FOR  EACH  MONTH. 

JAM.  -DEC.  > 

IH1ST0U2.  4.  97) "HISTOGRAM  VALUES 

THE  FIRST  SUBSCRIPT  VARIES  FROM  1  TO  12  AND  INDEXES  THE 
BAND  FOR  EACH  HISTOGRAM.  . 

THE  SECOND  SUBSCRIPT  INDEXES  THE  SEASON  OF  THAT  HISTOGRAM 

1  IMPLIES  WINTER  <DJF> 

2  IMPLIES  SPRING  (MAM) 

3  IMPLIES  SUMMER  (JJA) 

4  IMPLIES  FALL  (SON) 

THE  THIRD  SUBSCRIPT  DENOTES  WHICH  SQUARE  KM  AS  NOTED 
ABOVE 

COMMON/COLUMN/CR 10(970.  10. 2>.  NPTINC.  NR C INC.  NPTKM.  NRCKH,  XNTOT 
COMMON/RUFOUT/ IVC197) . NU(97>.  ZHAT( 12.  97).  IHISTO( 12.  4.  97) 

INTEGER *2  GRID 

DIMENSION  ICQUNT (32).  ZZER0O2.  32).  XMEAN ( 1 2 > •  NPT (  12).  BNDLMT(  13) 
PARAMETER  DUMMY- 1.  E-2 

DATA  BNDLMT/10.  E-5.  5.  E-5.  1.  E-4.  5.  E-4.  1.  E-3.  3.  E-3.  1.  E-2.  3.  E-2. 

1  1.  E-l,  5.  E-l.  1.  0,  5.  0.  10.  0/ 

DATA  (ZZEROd.  1).  1-1.  12>/12*DUHMY/ 

DATA  (ZZEROd.  2).  I-i,  12>/4*7.  3E-3.  2.  OE-2.3.  23E-2,  4*4.  3E-2. 

1  2.  6E-2.  7.  SE-3/ 

DATA  (ZZEROd.  3).  I-I.  l2)/4*7.  3E-3.  1.  13E-2.  1.  6E-2.  4*2  OE-2. 

1  1.  4E-2.7  SE-3/ 

DATA  (ZZEROd.  4).  I-I,  12)/4*1.  OE-2.  2.  SE-Z.  4.  OE-2.  4*3.  3E-2.  3.  23E-2. 
1  1 . OE-2/ 

DATA  (ZZEROd.  3).  1-1,  12>/12*1.  1/  - 

DATA  (ZZEROd.  6).  1-1.  12>/4*.  IS.  .  43. .  72,  4*1.  0,  .  39,  IS/ 

DATA  (ZZEROd.  7).  1-1,  12)/4*.  64.  .  78, .  91.  4*1.  03.  85.  .  64/ 

DATA  (ZZEROd,  8),  1-1,  12)/4*.  32.  .  40.  .  48.  4*.  S3.  .  44, .  32/ 

DATA  (ZZEROd.  9),  I-I,  I2>/4*2.  4E-2.  3.  8E-2,  9.  3E-2.  4*1.  31-1, 8  7E-2. 

1  2.  4E-2/ 
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6000 

6090 

6100 

6190 

6200 

6290 

6300 

6390 

6400 

6900 

6990 

6600 

6690 

6700 

6800 

6900 

7000 

7100 

7200 

7300 

7400 

7500 

7600 

7700 

7800 

7900 

8000 

8100 

8200 

8300 

8400 

8900 

8600 

8700 

8800 

8900  10 

9000 
9100 
9200 

9300  20 

9400 
9900 
9600 

9700  30 

9800 

9900 

10000  C 
10100 
10200 
10300 
10400 
10900  C 
10600 
10700 
10800  C 
10900 
11000 
11100  C 
11200 
11300 
11400 


SATA  (ZZEROtl.  10).  >1.  12>/4*1.  49E-2.  3  6E-2,  5  BE-2.  4*8  0E-2.  4  3E-2 
1  1.  45E-2/ 

SATA  (ZZEROtl.  ID.  >1.  12)/4*2.  9E-2.  6  3E-2.  10.  2E-2.4*!.  4E-1.8  2E-2 
1  8.  9E-8/ 

DATA  (ZZEROtl.  12).  I-l,  12>/4*1  6E-2.  4  OE-2. 6  4E-2.  4*8.  OE-2.  4.  8E-2, 
1  1.6E-2/ 

DATA  < 2 ZERO ( 1.  13).  1*1.  12)/4*  9E-2.  1  29E-2.  2  OE-2.  4*2.  79E-2. 1  6E-2 
1  '  .  5E-2/ 

DATA  (ZZEROtl.  14). 1*1. 12)/12*E-2/ 

DATA  <  ZZEROt 1.  15). 1-1.  12) /4*1.  8E-4.  4  9E-4.  7.  2E-4.  4*1.  OE-3.  9.  9E-4. 
1  1  BE -4/ 

DATA  <  ZZEROt  1.  16).  1-1.  12>/4*1.  7E-3.  4  3E-3.  6.  9E-3.  4*9.  9E-3.  9.  6E-3. 
1  1.7E-3/ 

DATA  ( ZZEROt 1.  17).  1-1.  12>/4*7  5E-3. 2.  19E-2.  3  6E-2. 

1  4*9.  OE-2.  2.  9E-2. 7  5E-3/ 

DATA  < ZZEROt 1.  18). 1-1.  12)/12*9.  OE-4/ 

DATA  (ZZEROtl.  19). I-l.  12)/12*DUMNY/ 

DATA  (ZZEROtl.  20). 1-1.  12>/12*DUHHV/ 

DATA  ( ZZEROt  1.21).  1-1. 12)/12*DUKMY/ 

DATA  (ZZEROtl. 22). 1-1. 12>/12*DUMHY/ 

DATA  (ZZEROtl.  23). 1-1. 12>/12*DUMMY/ 

DATA  (  ZZEROtl.  24).  1-1.  12)/12*DUMMY/ 

DATA  (ZZEROtl.  29),  I-l.  12>/12»DUHMY/ 

DATA  ( ZZEROt 1.24). I-l. 12>/12*DUMMV/ 

DATA  ( ZZEROt I.  27). 1-1.  12)/12*DUNMY/ 

DATA  (ZZEROtl.  28). I-l.  12>/12*4.  OE-1/ 

DATA  (ZZEROtl.  29). 1-1.  12>/12*9.  OE-4/ 

DATA  (ZZEROtl.  30). I-l.  12>/12*4.  OE-1/ 

DATA  (ZZEROtl.  31 ). I-l, 12)/12*5.  9E-1/ 

DATA  (ZZEROtl.  32). I-l. 12)/12*B.  9E-1/ 

DO  1000  J-1.97 
NUt J)-0 
IMAX— 0 

DO  10  1*1.32 
1C0UNT 1 1 )-0 
CONTINUE 
DO  20  L-1.4 
DO  20  K— 1.  12 
IHlSTOtK.  L.  J)-0 
CONTINUE 
DO  30  1-1.  12 
XMEANt I )— 3.  0 
NPTt I )-0 
CONTINUE 
DO  90  L-l.NRCHH 
DO  80  K-l.NRTKN 

EXTRACT  SURFACE  TYRE.  INCREMENT  PROPER  COUNTER 
ISURF-ORIDt  ( J-l >*NPTKM*K.  L.  2) 

INDEX- ISURF+1 

I COUNT ( INDEX >-ICOUNT( INDEX >*1 
DO  90  I-l.  12 

IF  UNDEFINED  CODE  DO  NOT  PROCESS 

IFt  I  SURF.  EO.  0.  OR.  ( ISURF.  OE.  18.  AND.  ISURF.  LE.  26)  >00  TO  90 
HElQHT-ZZEROtl.  INDEX) 

COMPUTE  AVERAGE  AND  F1U.  HISTOGRAM 
XMEANt I ) -XMEANt 1 )+LOO<  HEIGHT) 

»*TtI>-NPT(I)*l 
GET  SEASON  INDEX 
IFt  1.  EO.  12)  THEN 
ISEA8N— 1 

ELSE 
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1 1  900 
11600 
11700 

neoo 

11900 

12000 

12100 

12200 

12300 

12400 

12900 

12600 

12700 

12800 

12900 

13000 

13100 

13200 

13300 

13400 

13300 

13600 

13700 

13800 

13900 

14000 

14100 

14200 

14300 

14400 

14900 

14600 

14700 

14800 

14900 

19000 

13100 

19200 


1 


40 

49 

90 

80 

90 


100 


1SEABN-1/3+1 

END  IF 

GET  BAND  INDEX  * 

DO  40  H-I. 1$ 

IF1HEICHT  OE.  INDLNT(H)  AND  HE16HT.  LT  BNDLNTUI+l > )  THEN 
IBAND-H 
CO  T3  49 
END  IF 
CONTINUE 
STOP  'CETRUF ' 

1H1ST01 1 BAND-  ISEASN.  J)»IHISTO( I  BAND.  ISEA8N.  J)  +  l 
CONTINUE 
CONTINUE 
CONTINUE 
DO  100  I-l.  12 

XMEAN ( I ) ■ XNEAN ( I > /FLOAT < NFT ( I > ) 

2HATCI.  J)-EXP<XMEAN(X>> 

CONTINUE 
DO  110  1-1.32 
IF <  I.  EO.  1>  THEN 

NU( J)-NU<  J) +1 COUNT ( I > 

ELSE  IF<1.  OE  2.  AND.  I.  LE.  IS)  THEN 
XF< I COUNT < I ) .  OT.  I MAX)  THEN 
IVC< J>-1-1 
I  MAX- 1 COUNT! I) 

ENOIP 

ELSE  IFd.OE.  19.  AND.  I.  LE.  27)  THEN 
NU<J)—NU(U)+1 COUNT ( 1 ) 


ELSE 


110 

lOOO 


EN01F 

CONTINUE 

CONTINUE 

RETURN 

END 


I F ( 1 COUNT ( 1 ) .  CT.  I MAX) 
IVC(J)-1-1 
INAX— ICOUNTt I > 
END  IF 


THEN 


) 

1 
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ICO 

PROGRAM  TDRUFF 

200 

C-— — 

300 

C-  THIS 

FROGRAM  COMPUTES  VALUES  OF  TERRAIN-DERIVED  AERODYNAMIC  ROUGHNESS 

400 

C-  AND 

TERRAIN  STATISTICS  FOR  ELEVATION  DATA  SUPPLIED  ON  A  SQUARE  GR ID 

500 

C- 

600 

C-  TDRUFF  HAS  CHEATED  FOR  THE  US  ARMY  ATMOSPHERIC  SCIENCES  LABCRAT  -v  by 

700 

C-  THE 

UNIVERSITY  OF  DAYTON  RESEARCH  INSTITUTE  31  MARCH  1981 

300 

C- 

900 

C-  NOTE 

:  SOME  PARTS  OF  THIS  CODE  MAY  BE  SPECIFIC  TO  DEC  VAX/VMS  FCF'S-n 

i  000 

^  _ 

1 100 

C- 

1200 

INTEGERS  1DATA(2300.  25).  I WORK <23.  23) 

:300 

DATA  I DATA/62300  *  -1000/ 

1400 

DATA  I WORK/ 623 4-1 000/ 

1500 

c- 

1600 

C-  READ 

IN  SCAN  LINES  TO  FORM  A  COLUMN  OF  1HM  WIDE  DATA 

1700 

C- 

1800 

CPEN(UNIT-2,  RECORDTYPE* 'VARIABLE'.  STATUS- 'OLD ' , 

1900 

X 

FORM- 'UNFORMATTED'.  RECORDSI 2E-7300) 

2000 

C- 

2100 

C-  READ 

PAST  THE  TWO  HEADER  RECORDS  AT  THE  BEGINNING  CF  THE  TAPE 

2200 

C- 

2300 

READ (2) 

2400 

READ (2) 

2300 

I  EOF  -  0, 

2600 

IX  -  0 

2700 

I COL  -  1 

2800 

IREC  -  0 

2900 

READ  •.  BASEN.  BASSE.  NPPCOL. DELTA 

2000 

STARTE  -  FLOAT < 1PIX( BASSE)  /  1000)  *  lOOO. 

3100 

STARTN  -  FLOAT! IF IX (BASEN)/ 1000)  *  1000. 

2200 

ENDNOR  -  STARTN  *  1000. 

3300 

ENDEAS  -  STARTE  ♦  lOOO. 

3400 

ELAST  «  BASSE  -  DELTA 

3300 

ENEXT  -  ELAST  ♦  DELTA 

3600 

PRINT  *.  'BASEN.  BASEE.  NPPCOL.  DELTA:  ' ,  BASEN. BA5EE. NPPCOL.  DELTA 

3700 

C- 

2800 

C-  NPPCOL  -  *  OF  POINTS  PER  SCAN  LINE ( COLUMN ) 

3900 

C- 

4000 

50 

CONTINUE 

4100 

IF<ENEXT.  LT.  ENDEAS)  THEN 

4200 

READ <2.  END-32)  (IDATAdROW.  ICOL).  I ROW- 1.  NPPCOL) 

4300 

GO  TO  54 

4400 

52 

I EOF  -  1 

4300 

00  TO  33 

4600 

54 

CONTINUE 

4700 

ELAST  -  ENEXT 

4800 

ENEXT  -  ENEXT  «■  DELTA 

4900 

IREC  -  IREC  -  1 

3000 

ICOL  -  ICOL  a  1 

3100 

00  TO  50 

5200 

ENDIF 

3300 

C- 

3400 

C-  A  1  KM  HIDE  COLUMN  HAS  BEEN  READ  IN.  NOW  DIVIDE  IT  INTO 

3500 

C-  1KM  TALL  SECTORS  OHM  BY  1KM>  FOR  FURTHER  PROCESSING 

3600 

C- 

5700 

53 

ICOL  -  ICOL  -  1 

5800 

C- 

5900 

C-  I COL* 

0  OCCURS  WHEN  THE  PROORAM  ENCOUNTERS  AN  END  OF  FILE  WHEN  READIKO 

6000 

C-  THE 

FIRST  RECORD  OF  THE  NEXT  COLUMN  OF  SQUARE  KILOMETERS. 

6100 

C- 

35 


a  200 

6300 

6400 

1 900 

6600 

6700 

6800 

6900 

7000 

70 

7100 

73 

7200 

7300 

7400 

79 

7900 

7a00 

7700 

7800 

7<»00 

3000 

C- 

3100 

C- 

3200 

c- 

9300 

c- 

8400 

S900 

9600 

9700 

c- 

8800 

c- 

3900 

c- 

9000 

c- 

9100 

c- 

9200 

c- 

9300 

9400 

c- 

9900 

c- 

9600 

c- 

9700 

c- 

9800 

9900 

lOOOO 

10100 

10200 

10300 

c- 

10400 

c- 

10900. 

c- 

10600 

c- 

10700 

c- 

10800 

c- 

10900 

c- 

11000 

c- 

11100 

11200 

11300 

11400 

11900 

11600 

11700 

11800 

80 

11900 

12000 

c- 

12100 

c- 

12200 

c- 

1F< I COL. EO.  0)  60  TO  100 
IROWABS  •  1 

TOP  NOR  ■  FLOAT  <  I F I X  <  B A&EN+FLOAT  (  NPPCOL )  •DELTA  )  / 1 OOO  •  10C0>  * 

X  1000. 

STARTN  -  FLOAT < 1FIX ( BASEN) / 1000  *  1000) 

ENDNOR  »  STARTN  ♦  1000 

mow  -  i 

R OWL A ST  -  BASEN  -  DELTA 
ROWNEXT  -  ROWLAST  »  DELTA 
CONTINUE 

IF(ROWNEXT.  LT  ENDNOR)  THEN 
DO  79  1C  *  1,  I COL 

!WGRK<  mow.  1C)  •  IDATAC IROWABS.  IC> 

mow  -  mow  ♦  i 

IROWABS  ■  1R0WABS  ♦  1 
ROWLAST  -  ROWNEXT 
GO  TO  70 
ENDIF 

DONE  FILLING  THE  1W0RK  ARRAY < THIS  ARRAY  REPRESENTS 
A  1KH  BY  1KH  AREA).  ARRAY  HAS  DIMENSIONS  I ROW- 1.  I COL 

1R0W  -  1R0W  -  1 
IN  -  IFIX (STARTN) 

IE  •  IFIX (STARTE) 

SUBROUTINE  CALC  COMPUTES  THE  VARIOUS  VALUES  FOR  THE  SQUARE  KILOMETER 
REPRESENTED  IN  THE  ARRAY  IWORK. 

IN  IS  THE  NORTHING  OF  THE  LOWER  LEFT  CORNER  OF  THE  SQUARE 
IE  IS  THE  EAST  I  NO  OF  THE  LOWER  LEFT  CORNER  OF  THE  SQUARE. 

CALL  CALC < IWORK. 1C0L.  1R0W.  IN.  IE) 

SUBROUTINE  COMDAT  SIMPLY  RF-INITIALIZES  THE  IWORK  ARRAY  TO  NO-DATA 
(-1000  VALUES). 

CALL  COMDAT < IWORK. IROW. ICOL.  STARTN.  STARTE) 

STARTN  -  FLOAT ( IFIX (ENDNOR+l )/ 1000*1000) 

ENDNOR  -  STARTN  ♦  10C0 
IROW  -  1 

IF < STARTN.  LT.  TOPNOR)  00  TO  73 

THIS  ENDS  THE  DATA  COLUMN  PROCESSING  SECTION 
COMPUTE  PARAMETERS  TO  READ  IN  ANOTHER  1KH  WIDE 
COLUMN  OF  DATA 


-  TEST  FOR  END  OF  FILE  • 

• 

IFdEDF  EO.  1)  GO  TO  100 
STARTE  -  ENDEAS 
ENDEA8  -  ENDEAS  4-1001 

ENDEAS  -  FLOAT(dFIX(ENDEAS>d  1/1000*1000) 
ICOL  -  I 
DO  SO  1-1.2900 
DO  SO  J  -  1.  29 

IDATA(I.J)  -  -1000 

GO  TO  90 

END  OF  FILE  SECTION 
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pj  W  M  M 


:3oo 

too 

WR1TE<6> 110)  IREC 

!400 

110 

FORMAT (  '  RECORDS  READ  -  MlO) 

1900 

STOP  'END  OF  JOB' 

600 

END 

lOO 

SUBROUTINE  COMDAT < IWORK.  1RD1M.  ICDIM,  STARTN,  STARTE) 

300 

300 

INTECER*2  X WORK ( 39.  39 1 

400 

c- 

500 

c- 

RE~ INITIALIZE  THE  I WORK  ARRAY  TO  NO-DATA  VALUES 

600 

c- 

AREA  REPRESENTED  IN  THE  I WORK  ARRAY. 

700 

.  c- 

800 

DO  30  I  -  1.39 

900 

DO  30  J  ■  1.39 

1000 

30 

I WORM  I,  J>  -  -1000 

1100 

RETURN 

1200 

END 

100 

200 

200 

400 

500 

600 

700 

aoo 

900 

1000 

1100 

1200 

1300 

1400 

1500 

1600 

1700 

1800 

1900 

2000 

2100 

2200 

2300 

2400 

2500 

2600 

2700 

2800 

2900 

3000 

3100 

3200 

3300 

3400 

3500 

3600 

3700 

3800 

3900 

4000 

4100 

4200 

4300 

4400 

4500 

4600 

4700 

4800 

4900 

5000 

5100 

5200 

5300 

5400 

5500 

5600 

5700 

5800 

5900 

6000 

6100 


c- 


c- 


V- 

c- 

c- 

c- 

c- 


c- 

c- 

c- 

c- 

c- 

£- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 
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SUBROUTINE  CALC ( I DATA.  IXDIH.  IYDIM,  NORTH.  1 EAST ) 


this  subroutine  COMPUTES  THE  OUTPUT  VARIABLES  from  the  elevation 
DATA  FOR  EACH  OF  THE  SQUARE  KILOMETERS  CONTAINED  IN  AN  ARRAY  GRID 
IT  ALSO  CALLS  ROUTINE  PEAKS  TO  COUNT  THE  NUMBER  OF  PEAKS  AND  VALLEY 
IN  A  SCAN  LINE  ACROSS  1,  KM  SOU  ARE 

INPUT: 


IXDIH  -  THE  NUMBER  OF  POINTS  IN  THE  N-S  DIRECTION  PER  KM. 

IYDIM  ■  THE  NUMBER  OF  POINTS  IN  THE  E-W  DIRECTION  PER  KM. 

NORTH  -  NORTHING  OF  S-W  CORNER  POINT 

I  EAST  a  EASTING  OF  S-W  CORNER  POINT 

I DATA!  >  -  THE  FIRST  SUBSCRIPT  VARIES  FROM  1  TO  16  AND  INDEXES  THE 
POINTS  OF  1  SQUARE  KM  IN  THE  N-S  DIRECTION. 

THE  SECOND  SUBSCRIPT  VARIES  FROM  1  TO  16  AND  INDEXES  THE 
POINTS  OF  1  SQUARE  KM  IN  THE  E-W  DIRECTION. 

OUTPUT: 


NORTH  >  NORTHING  OF  S-U  CORNER  POINT. 

IE AST  a  EASTING  OF  S-U  CORNER  POINT. 

HBAR  -  THE  MEAN  ELEVATION  (METERS). 

SIGMA  -  THE  ELEVATION  STANDARD  DEVIATION  (METERS) 

HOBAR  -  THE  TERRAIN  ROUGHNESS  PARAMETER  (METERS). 
ICOUNT  -  NUMBER  OF  ELEVATION  DATA  POINTS  IN  A  SQUARE  KM. 
IMIN  -  THE  MINIMUM  ELEVATION  (METERS) 

I MAX  -  THE  MAXIMUM  ELEVATION  (METERS). 


INTECER*2  IDATA(2S.  25).  IARAY(30> 

SUM  »  O. 

I SCANS  a  o 
SUMSQ  a  0. 

I MAX  a  O 
IMIN  -  30000 
NPV  a  0 
HOBAR  a  O.  - 
ICOUNT  -  0 

DO  SO  IROW  -  1. IYDIM 
DO  SO  I COL  a  l,  IXDIM 

IF( I DATA (IROW.  ICOL).  EO.  -lOOO)  00  TO  SO 

ICOUNT  -  ICOUNT  ♦  1 

DATA  -  FLOAT (I DATA (I ROW. ICOL) > 

SUM  ■  SUM  ♦  DATA 

SUMSQ  -  SUMSQ  a  ( DATA*OATA) 

IF< I DATA ( IROW.  ICOL).  LT.  IMIN)  IMIN  -  I DATA! IROW.  ICOL) 
IF( IDATA< IROW. ICOL).  OT.  IMAX)  IMAX.-  I DATA ( IROW. ICOL) 
CONTINUE 

IF< ICOUNT.  EQ.  0)  RETURN 
HBAR  -  SUM  /  FLOAT < ICOUNT ) 

IF( ICOUNT.  EQ.  1)  THEN 
SIGMA  -  0. 

HOBAR  -  0. 

GO  TO  300 
END  IF 

A  -  SUMSQ  -  (3UM*?UM/FL0AT( ICOUNT) ) 

IF (  A.  LT.  0.  )  THEN 
SIGMA  -  o. 

ELSE 
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SIGMA  -  SORT < A  /  FLOAT! I COUNT-1  > > 


o200 

6300 

6400 

6300 

s600 

6700 

6300 

a<?00 

7000 

7100 

7200 

7300 

7400 

7300 

7600 

7700 

7800 

7900 

3000 

3100 

3200 

3300 

8400 

9300 

8600 

8700 

8300 

3900 

9000 

9100 

9200 

9200 

9400 

9300 

9600 

9700 

9800 

9900 

lOOOO 

10100 

10200 

10300 

10400 

10300 

10600 

10700 

10800 

10900 

11000 


ENDIF 

c- 

C-  FOR  EACH  ROW  .BUILD  AN  ARRAY  OF  DATA  VALUES 
C- 

DO  100  IRC W  -  1.  IYDIH 
IC  -  0 

DO  90  1C0L  -  1. IXDIM 

IF ( I DATA! I ROW,  ICOL >  EO.  -1000)  CO  TO  90 
IC  -  IC  ♦  1 

lARAY(IC)  -  IDATAI IRCW. ICOL) 

90  CONTINUE 

C- 

C-  COMPUTE  THE  NUMBER  OF  PEAKS  AND  VALLEYS  FOR  THIS  ROW 
C- 

CALL  PEAKS < IARAY, IC,  HO.  NPNV,  1SCANS.  ISCANS) 

HOBAR  -  HOBAR  «■  HO 
100  NPV  “  NPV  «■  NPNV 

c- 

C-  FOR  EACH  COLUMN  BUILD  AN  ARRAY  OF  DATA  VALUES 
C- 

DO  200  ICOL  -  1, IXDIM 
IC  -  0 

DO  190  1RQW  »  1, IYDIM 

IF< IDATAI IROW. ICOL).  EO. -1000)  GOTO  190 
IC  •  IC  +  1 

IARAY ( 1C )  -  IDATAI IROU.  ICOL) 

190  CONTINUE 

C- 

C-  COMPUTE  THE  NUMBER  OF  PEAKS  AND  VALLEYS  FOR  THIS  COLUMN 
C- 

CALL  PEAKS! IARAY.  IC.  HO,  NPNV.  ISCANS) 

HOBAR  -  ilOBAR  ♦  HO 
200  NPV  *  NPV  ♦  NPNV' 

C-  . 

C-  COMPUTE  FINAL  VALUES  AND  OUTPUT  THEM 
C- 

300  IF! ISCANS.  EQ.  0)  THEN 

HOBAR  *  0 

ELSE 

HOBAR  -  HOBAR  /  FLOAT! ISCANS) 

ENDIF 

C- 

C-  WRITE  THE  OUTPUT  RECORD  FOR  THIS  KILOMETER  SQUARE. 

•  C-' 

WRITE (7.  1000)  NORTH,  I EAST .  HBAR,  SIGMA.  HOBAR.  ICOUNT.  IMIN.  IMAX 
1000  F0RMATI2IS.  3E14.  7. 318) 

RETURN 

END 


SUBROUTINE  PEAKS ( IARAY, IC,  HO.  NPNV, ISCANS) 

INTEGER*:  I ARAY < 30) 

THIS  ROUTINE  DETERMINES  THE  NUMBER  OF  PEAKS  AND  VALLEYS  IN  A  SCAN 
LINE  ACROSS  A  1  KM  SQUARE  AND  THE  TDR  (HO)  FOR  THE  SCAN  LINE 
(TRAVERSE).  THIS  SCAN  LINE  IS  STORED  IN  I ARAY. 

NPNV  «  0 
NPEAKS  =»  0 
NVALLS  -  0 

HERE  is  THE  DEFAULT  SEiTING  FOR  HO 

HO  »  1.  c-2 
IF ( IC.  L£.  2)  RETURN 
ISCANS  =  ISCANS  *  1 
DTOTAL  -  0. 

IF( I ARA * ( 1 ) .  CT.  I ARAY ( 2) )  THEN 
NPEAKS  -  NPEAKS  ♦  1 

ELSE 

IF( IARAY ( 1 )  LI .  IARAYO)  THEN 
NVALLS  *  NVALLS  ♦  t 

END  IF 

END  IF 

ICM1  -  IC  -  1 
DO  30  I  -  2.  ICM1 

IF( ( IARAY ( I ) .  GT  lARAY(I-l))  AND.  < I ARAY ( I > .  GT. 

X  IARAY ( 1+1 ) ) )  THEN 

NPEAKS  *  NPEAKS  *  1 

ELSE 

IF ( ( IARAY  < I ) .  LT.  IARAY (I  —  X  >  > .  AND.  ( IARAY ( I > . LT. 

X  IARAY ( 1*1 >) )  THEN  ' 

NVALLS  *  NVALLS  ♦  1 

ENDIF 

ENDIF 

DTOTAL  -  DTOTAL  ABS( FLOAT (  IARAY  ( I ) -IARAY  ( 1-1 )  )  ) 


IF( I ARAY ( IC).  GT.  IARAY(IC-l))  THEN 
NPEAKS  -  NPEAKS  *  1 

ENDIF 

IF( I ARAY ( IC).  LT.  IARAY (IC-!))  THEN 
NVALLS  -  NVALLS  ♦  1 

ENDIF 

DTOTAL  •  DTOTAL  ♦  ABS ( FLOAT < IARAY ( IC >  -  IARAY ( IC-1 ) ) ) 
NPNV  *  NPEAKS  ♦  NVALLS 
IF (NPNV.  LE.  1)  THEN 
DTOTAL  -  0. 

HO  -  I.  E-2 

ELSE 

HO  -  (DTOTAL  *  DTOTAL)  /  (8000.  *  (NPNV-D) 

ENDIF 

RETURN 

EM) 


